# Extreme weather events do not increase political parties' environmental attention

#'   This scirpt
#'    > performs PanelMatch models based on relief classification
#'    > creates one graph for all models
#'    > exports Extended Data Figure 5



# Tim Wappenhans, António Valentim, Heike Klüver, and Lukas F. Stoetzer
# First: 22-02-2023
# Last: 10-03-2024


# PACKAGES --------------------------------------------------------------------
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  tidyverse,
  PanelMatch,
  panelView,
  ggpubr,
  ggrepel,
  cowplot,
  here
)


# DATA WRANGLING FOR PANELMATCH ------------------------------------------------
df.panelmatch <- read_csv(here("data/data_output/press_weekly.csv"))   |> 
  mutate(
    
    # only integers allowed 
    parlgov_id = as.integer(parlgov_id),
    t = as.integer(t),
    country = unclass(as.factor(country_name)),
    family = unclass(as.factor(family_name_short))
    
  ) |> 
  
  select(
    
    # Covs
    parlgov_id, t, family, country, cabinet_party,
    
    # Treatments
    treat, treat_flood, treat_storm, treat_fire, treat_temp,
    
    # Outcome
    climate_prs, relief_prs, relief_words, pr_total
    
  ) |> 
  
  mutate(
    
    # relief
    outcome_relief = relief_prs / pr_total,
    outcome_relief = ifelse(is.nan(outcome_relief), 0, outcome_relief),

    outcome_relief_words = relief_words,
    outcome_relief_words = ifelse(is.nan(outcome_relief_words), 0, outcome_relief_words)
    
  ) |> 
  
  as.data.frame() |> 
  
  filter(!is.na(parlgov_id)) 


# Globally define leads and lags, iterations
lag <- 6
lead <- 6
n_iterations <- 1000


# create empty tibble for data viz
gg.did <- data.frame(
  coef = numeric(), 
  se = numeric(),
  model = numeric()
)


# LOOPS TO CREATE ONE GRAPH ACROSS ALL MODELS ----------------------------------
models <- c("Relief PRs",
            "Relief Words")

for (model in models) {
  if (model == "Relief PRs") {
    
    PM.results <- PanelMatch(lag = lag, 
                             lead = 0:lead, 
                             data = df.panelmatch,
                             time.id = "t", 
                             unit.id = "parlgov_id",
                             treatment = "treat", 
                             outcome.var = "outcome_relief",
                             match.missing = TRUE,
                             refinement.method = "mahalanobis",
                             covs.formula = ~ I(lag(outcome_relief, 1:lead)),
                             exact.match.variables = "family",
                             size.match = 10, 
                             qoi = "att",
                             forbid.treatment.reversal = FALSE)
    
  } else if (model == "Relief Words") {
    
    PM.results <- PanelMatch(lag = lag, 
                             lead = 0:lead, 
                             data = df.panelmatch,
                             time.id = "t", 
                             unit.id = "parlgov_id",
                             treatment = "treat", 
                             outcome.var = "outcome_relief_words",
                             match.missing = TRUE,
                             refinement.method = "mahalanobis",
                             covs.formula = ~ I(lag(outcome_relief_words, 1:lead)),
                             exact.match.variables = "family",
                             size.match = 10, 
                             qoi = "att",
                             forbid.treatment.reversal = FALSE)
  }
  
  
  ## Treatment
  PE.treat <- PanelEstimate(sets = PM.results, 
                            data = df.panelmatch,
                            se.method = "bootstrap",
                            number.iterations = n_iterations)
  
  # Store Coefs and SEs
  gg.treat <- summary(PE.treat, 
                      verbose = F, 
                      bias.corrected = F) |> 
    as.data.frame() |> 
    mutate(model = model,
           period = row_number() - 1) |> 
    rename(low_ci = `2.5%`,
           up_ci = `97.5%`) |> 
    select(-std.error)
  
  ## Placebo
  PE.placebo <- placebo_test(PM.results, 
                             data = df.panelmatch, 
                             number.iterations = n_iterations)
  
  # placebo estimates  
  placebo_estimates <- PE.placebo$estimates |> as_tibble() |> 
    mutate(period = row_number() - (lead + 1)) |> 
    rename(estimate = value) 
  
  
  # placebo CIs, percentiles from bootstrapping
  percentiles <- PE.placebo$bootstrapped.estimates |> as_tibble() |> 
    summarize(across(everything(), ~ quantile(., probs = c(0.025, 0.975))))
  
  percentiles_longer <- percentiles |> 
    mutate(quantile = row_number(),
           quantile = case_when(quantile == 1 ~ "low_ci",
                                quantile == 2 ~ "up_ci")) |> 
    pivot_longer(cols = starts_with("t-"), names_to = "period", values_to = "value") %>%
    pivot_wider(names_from = "quantile", values_from = "value") |> 
    mutate(period = substr(period,2,3),
           period = as.numeric(period))
  
  # join placebo estimates and placebo CIs
  gg.placebo <- left_join(placebo_estimates, percentiles_longer,
                          by = "period") |> 
    mutate(model = model)
  
  # bind placebo to treatment
  gg.tmp <- rbind(gg.placebo, gg.treat)
  
  # Add to main viz DF
  gg.did <- rbind(gg.did, gg.tmp)
  
}


# PLOTS ------------------------------------------------------------------------
# globally define settings
pre_color <- "#4885C1"
post_color <- "#AE3A4E"

pre_fill <- "#4885C1"
post_fill <- "#AE3A4E"


pre_shape <- 17
post_shape <- 16

# create indicator
gg.did <- gg.did |> 
  mutate(pre = ifelse(period < 0, "pre", "post"))

# reorder
gg.did$model <- factor(gg.did$model, levels = models) 



# Store it so don't have to bootstrap every time anew
write_csv(gg.did,
          file = here("data/data_output/ggdid_dic_relief.csv"))

# pick back up
gg.did <- read_csv(here("data/data_output/ggdid_dic_relief.csv"))

# data wrangling for plot
gg.did <- gg.did |> 
  mutate(
    
    # make it percentage points
    estimate = estimate * 100, 
    low_ci = low_ci * 100, 
    up_ci = up_ci * 100
    )



gg.did |> 
  ggplot() +
  
  # Gray pre treatment shade
  geom_rect(aes(xmin = -Inf, xmax = -0.5, ymin = -Inf, ymax = Inf),
            fill = "#e0e0e0", alpha = 0.2) +
  
  # CIs
  geom_linerange(aes(x = period, ymin = low_ci, ymax = up_ci, color = pre),
                 alpha = .7,
                 linewidth = 0.75) + 
  
  # Coef
  geom_point(aes(x = period, y = estimate, shape = pre, color = pre),
             size = 3) + 
  
  # fixed lines
  geom_hline(yintercept = 0, linetype = 'dotted') + 
  
  # look
  ylab("") +
  xlab("Weeks relative to extreme weather event") + 
  scale_x_continuous(breaks = gg.did$period, labels = gg.did$period) +
  facet_wrap(~ model, 
             scales = "free", 
             ncol = 2) +
  scale_color_manual(values = c("pre" = pre_color,
                                "post" = post_color)) +
  scale_shape_manual(values = c("pre" = pre_shape, 
                                "post" = post_shape)) +
  scale_fill_manual(values = c("pre" = pre_fill, 
                               "post" = post_fill)) +
  guides(shape = "none",
         fill = "none",
         color = "none") +
  
  # theme
  theme_bw() +
  theme(panel.grid = element_blank(),
        plot.title = element_text(hjust = 0.5)) 

  
# Export
ggsave(here("results/graphs/fig_ex_data_05.pdf"),
       width = 10,
       height = 5)